home *** CD-ROM | disk | FTP | other *** search
/ Power Programmierung 2 / Power-Programmierung CD 2 (Tewi)(1994).iso / gnu / gnulib / sipp / libsipp / noise.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-07-28  |  9.0 KB  |  306 lines
  1. #include <math.h>
  2.  
  3. #include <sipp.h>
  4. #include <noise.h>
  5.  
  6. /* ================================================================ */
  7. /*                      Noise and friends                           */
  8.  
  9.  
  10. /*
  11.  *    Noise and Dnoise routines
  12.  *
  13.  *    Many thanks to Jon Buller of Colorado State University for these
  14.  *    routines.
  15.  */
  16.  
  17.  
  18. #define NUMPTS    512
  19. #define P1    173
  20. #define P2    263
  21. #define P3    337
  22. #define phi    0.6180339
  23.  
  24.  
  25. bool noise_ready = FALSE;
  26.  
  27. static double pts[NUMPTS];
  28.  
  29.  
  30. /*
  31.  * Doubles might have values that is too large to be
  32.  * stored in an int. If this is the case we "fold" the
  33.  * value about MAXINT or MININT until it is of an
  34.  * acceptable size.
  35.  *
  36.  * NOTE: 32 bit integers is assumed.
  37.  */
  38. static double
  39. iadjust(f)
  40.     double f;
  41. {
  42.     while ((f > 2147483647.0) || (f < -2147483648.0)) {
  43.         if (f > 2147483647.0) {         /* 2**31 - 1 */
  44.             f = (4294967294.0 - f);     /* 2**32 - 2 */
  45.         } else if (f < -2147483648.0){
  46.             f = (-4294967296.0 - f);
  47.         }
  48.     }
  49.     return f;
  50. }
  51.  
  52.  
  53.  
  54. /*
  55.  * Initialize the array of random numbers.
  56.  * RANDOM() is defined in sipp.h
  57.  */
  58. void noise_init()
  59. {
  60.     int i;
  61.    
  62.     for (i = 0; i < NUMPTS; ++i)
  63.         pts[i] = RANDOM();
  64.     noise_ready = TRUE;
  65. }
  66.  
  67.  
  68.  
  69. double noise(v)
  70. Vector *v;
  71. {
  72.    Vector p;
  73.    int xi, yi, zi;
  74.    int xa, xb, xc, ya, yb, yc, za, zb, zc;
  75.    double xf, yf, zf;
  76.    double x2, x1, x0, y2, y1, y0, z2, z1, z0;
  77.    double p000, p100, p200, p010, p110, p210, p020, p120, p220;
  78.    double p001, p101, p201, p011, p111, p211, p021, p121, p221;
  79.    double p002, p102, p202, p012, p112, p212, p022, p122, p222;
  80.    double tmp;
  81.  
  82.    p = *v;
  83.  
  84.    p.x = iadjust(p.x);
  85.    p.y = iadjust(p.y);
  86.    p.z = iadjust(p.z);
  87.    xi = floor(p.x);
  88.    xa = floor(P1 * (xi * phi - floor(xi * phi)));
  89.    xb = floor(P1 * ((xi + 1) * phi - floor((xi + 1) * phi)));
  90.    xc = floor(P1 * ((xi + 2) * phi - floor((xi + 2) * phi)));
  91.  
  92.    yi = floor(p.y);
  93.    ya = floor(P2 * (yi * phi - floor(yi * phi)));
  94.    yb = floor(P2 * ((yi + 1) * phi - floor((yi + 1) * phi)));
  95.    yc = floor(P2 * ((yi + 2) * phi - floor((yi + 2) * phi)));
  96.  
  97.    zi = floor(p.z);
  98.    za = floor(P3 * (zi * phi - floor(zi * phi)));
  99.    zb = floor(P3 * ((zi + 1) * phi - floor((zi + 1) * phi)));
  100.    zc = floor(P3 * ((zi + 2) * phi - floor((zi + 2) * phi)));
  101.  
  102.    p000 = pts[(xa + ya + za) % NUMPTS];
  103.    p100 = pts[(xb + ya + za) % NUMPTS];
  104.    p200 = pts[(xc + ya + za) % NUMPTS];
  105.    p010 = pts[(xa + yb + za) % NUMPTS];
  106.    p110 = pts[(xb + yb + za) % NUMPTS];
  107.    p210 = pts[(xc + yb + za) % NUMPTS];
  108.    p020 = pts[(xa + yc + za) % NUMPTS];
  109.    p120 = pts[(xb + yc + za) % NUMPTS];
  110.    p220 = pts[(xc + yc + za) % NUMPTS];
  111.    p001 = pts[(xa + ya + zb) % NUMPTS];
  112.    p101 = pts[(xb + ya + zb) % NUMPTS];
  113.    p201 = pts[(xc + ya + zb) % NUMPTS];
  114.    p011 = pts[(xa + yb + zb) % NUMPTS];
  115.    p111 = pts[(xb + yb + zb) % NUMPTS];
  116.    p211 = pts[(xc + yb + zb) % NUMPTS];
  117.    p021 = pts[(xa + yc + zb) % NUMPTS];
  118.    p121 = pts[(xb + yc + zb) % NUMPTS];
  119.    p221 = pts[(xc + yc + zb) % NUMPTS];
  120.    p002 = pts[(xa + ya + zc) % NUMPTS];
  121.    p102 = pts[(xb + ya + zc) % NUMPTS];
  122.    p202 = pts[(xc + ya + zc) % NUMPTS];
  123.    p012 = pts[(xa + yb + zc) % NUMPTS];
  124.    p112 = pts[(xb + yb + zc) % NUMPTS];
  125.    p212 = pts[(xc + yb + zc) % NUMPTS];
  126.    p022 = pts[(xa + yc + zc) % NUMPTS];
  127.    p122 = pts[(xb + yc + zc) % NUMPTS];
  128.    p222 = pts[(xc + yc + zc) % NUMPTS];
  129.  
  130.    xf = p.x - xi;
  131.    x1 = xf * xf;
  132.    x2 = 0.5 * x1;
  133.    x1 = 0.5 + xf - x1;
  134.    x0 = 0.5 - xf + x2;
  135.  
  136.    yf = p.y - yi;
  137.    y1 = yf * yf;
  138.    y2 = 0.5 * y1;
  139.    y1 = 0.5 + yf - y1;
  140.    y0 = 0.5 - yf + y2;
  141.  
  142.    zf = p.z - zi;
  143.    z1 = zf * zf;
  144.    z2 = 0.5 * z1;
  145.    z1 = 0.5 + zf - z1;
  146.    z0 = 0.5 - zf + z2;
  147.  
  148.    /*
  149.     * This expression is split up since some compilers
  150.     * chokes on expressions of this size.
  151.     */
  152.    tmp =  z0 * (y0 * (x0 * p000 + x1 * p100 + x2 * p200) +
  153.                 y1 * (x0 * p010 + x1 * p110 + x2 * p210) +
  154.                 y2 * (x0 * p020 + x1 * p120 + x2 * p220));
  155.    tmp += z1 * (y0 * (x0 * p001 + x1 * p101 + x2 * p201) +
  156.                 y1 * (x0 * p011 + x1 * p111 + x2 * p211) +
  157.                 y2 * (x0 * p021 + x1 * p121 + x2 * p221));
  158.    tmp += z2 * (y0 * (x0 * p002 + x1 * p102 + x2 * p202) +
  159.                 y1 * (x0 * p012 + x1 * p112 + x2 * p212) +
  160.                 y2 * (x0 * p022 + x1 * p122 + x2 * p222));
  161.  
  162.    return tmp;
  163. }
  164.  
  165.  
  166.  
  167. Vector Dnoise(p)
  168. Vector *p;
  169. {
  170.    Vector v;
  171.    int xi, yi, zi;
  172.    int xa, xb, xc, ya, yb, yc, za, zb, zc;
  173.    double xf, yf, zf;
  174.    double x2, x1, x0, y2, y1, y0, z2, z1, z0;
  175.    double xd2, xd1, xd0, yd2, yd1, yd0, zd2, zd1, zd0;
  176.    double p000, p100, p200, p010, p110, p210, p020, p120, p220;
  177.    double p001, p101, p201, p011, p111, p211, p021, p121, p221;
  178.    double p002, p102, p202, p012, p112, p212, p022, p122, p222;
  179.  
  180.    xi = floor(p->x);
  181.    xa = floor(P1 * (xi * phi - floor(xi * phi)));
  182.    xb = floor(P1 * ((xi + 1) * phi - floor((xi + 1) * phi)));
  183.    xc = floor(P1 * ((xi + 2) * phi - floor((xi + 2) * phi)));
  184.  
  185.    yi = floor(p->y);
  186.    ya = floor(P2 * (yi * phi - floor(yi * phi)));
  187.    yb = floor(P2 * ((yi + 1) * phi - floor((yi + 1) * phi)));
  188.    yc = floor(P2 * ((yi + 2) * phi - floor((yi + 2) * phi)));
  189.  
  190.    zi = floor(p->z);
  191.    za = floor(P3 * (zi * phi - floor(zi * phi)));
  192.    zb = floor(P3 * ((zi + 1) * phi - floor((zi + 1) * phi)));
  193.    zc = floor(P3 * ((zi + 2) * phi - floor((zi + 2) * phi)));
  194.  
  195.    p000 = pts[(xa + ya + za) % NUMPTS];
  196.    p100 = pts[(xb + ya + za) % NUMPTS];
  197.    p200 = pts[(xc + ya + za) % NUMPTS];
  198.    p010 = pts[(xa + yb + za) % NUMPTS];
  199.    p110 = pts[(xb + yb + za) % NUMPTS];
  200.    p210 = pts[(xc + yb + za) % NUMPTS];
  201.    p020 = pts[(xa + yc + za) % NUMPTS];
  202.    p120 = pts[(xb + yc + za) % NUMPTS];
  203.    p220 = pts[(xc + yc + za) % NUMPTS];
  204.    p001 = pts[(xa + ya + zb) % NUMPTS];
  205.    p101 = pts[(xb + ya + zb) % NUMPTS];
  206.    p201 = pts[(xc + ya + zb) % NUMPTS];
  207.    p011 = pts[(xa + yb + zb) % NUMPTS];
  208.    p111 = pts[(xb + yb + zb) % NUMPTS];
  209.    p211 = pts[(xc + yb + zb) % NUMPTS];
  210.    p021 = pts[(xa + yc + zb) % NUMPTS];
  211.    p121 = pts[(xb + yc + zb) % NUMPTS];
  212.    p221 = pts[(xc + yc + zb) % NUMPTS];
  213.    p002 = pts[(xa + ya + zc) % NUMPTS];
  214.    p102 = pts[(xb + ya + zc) % NUMPTS];
  215.    p202 = pts[(xc + ya + zc) % NUMPTS];
  216.    p012 = pts[(xa + yb + zc) % NUMPTS];
  217.    p112 = pts[(xb + yb + zc) % NUMPTS];
  218.    p212 = pts[(xc + yb + zc) % NUMPTS];
  219.    p022 = pts[(xa + yc + zc) % NUMPTS];
  220.    p122 = pts[(xb + yc + zc) % NUMPTS];
  221.    p222 = pts[(xc + yc + zc) % NUMPTS];
  222.  
  223.    xf = p->x - xi;
  224.    x1 = xf * xf;
  225.    x2 = 0.5 * x1;
  226.    x1 = 0.5 + xf - x1;
  227.    x0 = 0.5 - xf + x2;
  228.    xd2 = xf;
  229.    xd1 = 1.0 - xf - xf;
  230.    xd0 = xf - 1.0;
  231.  
  232.    yf = p->y - yi;
  233.    y1 = yf * yf;
  234.    y2 = 0.5 * y1;
  235.    y1 = 0.5 + yf - y1;
  236.    y0 = 0.5 - yf + y2;
  237.    yd2 = yf;
  238.    yd1 = 1.0 - yf - yf;
  239.    yd0 = yf - 1.0;
  240.  
  241.    zf = p->z - zi;
  242.    z1 = zf * zf;
  243.    z2 = 0.5 * z1;
  244.    z1 = 0.5 + zf - z1;
  245.    z0 = 0.5 - zf + z2;
  246.    zd2 = zf;
  247.    zd1 = 1.0 - zf - zf;
  248.    zd0 = zf - 1.0;
  249.  
  250.    /*
  251.     * This expressions are split up since some compilers
  252.     * chokes on expressions of this size.
  253.     */
  254.    v.x        = z0 * (y0 * (xd0 * p000 + xd1 * p100 + xd2 * p200) +
  255.                       y1 * (xd0 * p010 + xd1 * p110 + xd2 * p210) +
  256.                       y2 * (xd0 * p020 + xd1 * p120 + xd2 * p220));
  257.    v.x       += z1 * (y0 * (xd0 * p001 + xd1 * p101 + xd2 * p201) +
  258.                       y1 * (xd0 * p011 + xd1 * p111 + xd2 * p211) +
  259.                       y2 * (xd0 * p021 + xd1 * p121 + xd2 * p221));
  260.    v.x       += z2 * (y0 * (xd0 * p002 + xd1 * p102 + xd2 * p202) +
  261.                       y1 * (xd0 * p012 + xd1 * p112 + xd2 * p212) +
  262.                       y2 * (xd0 * p022 + xd1 * p122 + xd2 * p222));
  263.                                   
  264.    v.y        = z0 * (yd0 * (x0 * p000 + x1 * p100 + x2 * p200) +
  265.                       yd1 * (x0 * p010 + x1 * p110 + x2 * p210) +
  266.                       yd2 * (x0 * p020 + x1 * p120 + x2 * p220));
  267.    v.y       += z1 * (yd0 * (x0 * p001 + x1 * p101 + x2 * p201) +
  268.                       yd1 * (x0 * p011 + x1 * p111 + x2 * p211) +
  269.                       yd2 * (x0 * p021 + x1 * p121 + x2 * p221));
  270.    v.y       += z2 * (yd0 * (x0 * p002 + x1 * p102 + x2 * p202) +
  271.                       yd1 * (x0 * p012 + x1 * p112 + x2 * p212) +
  272.                       yd2 * (x0 * p022 + x1 * p122 + x2 * p222));
  273.                                   
  274.    v.z        = zd0 * (y0 * (x0 * p000 + x1 * p100 + x2 * p200) +
  275.                        y1 * (x0 * p010 + x1 * p110 + x2 * p210) +
  276.                        y2 * (x0 * p020 + x1 * p120 + x2 * p220));
  277.    v.z       += zd1 * (y0 * (x0 * p001 + x1 * p101 + x2 * p201) +
  278.                        y1 * (x0 * p011 + x1 * p111 + x2 * p211) +
  279.                        y2 * (x0 * p021 + x1 * p121 + x2 * p221));
  280.    v.z       += zd2 * (y0 * (x0 * p002 + x1 * p102 + x2 * p202) +
  281.                        y1 * (x0 * p012 + x1 * p112 + x2 * p212) +
  282.                        y2 * (x0 * p022 + x1 * p122 + x2 * p222));
  283.  
  284.    return v;
  285. }
  286.  
  287.  
  288.  
  289. double turbulence(p, octaves)
  290.     Vector *p;
  291.     int     octaves;
  292. {
  293.     Vector tmp;
  294.     double scale = 1.0;
  295.     double t     = 0.0;
  296.  
  297.     while (octaves--) {
  298.         tmp.x = p->x * scale;
  299.         tmp.y = p->y * scale;
  300.         tmp.z = p->z * scale;
  301.         t += (noise(&tmp) / scale);
  302.         scale *= 2.0;
  303.     }
  304.     return t;
  305. }
  306.